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SUMMARY 

This report describes a new algorithm for fast, automatic Integration of 
chemical kinetic rate equations describing homogeneous, gas-phase combustion 
at constant pressure. Particular attention Is paid to the distinguishing 
physical and computational characteristics of the Induction, heat-release and 
equilibration regimes. The two-part predictor-corrector algorithm, based on 
an "exponentially-fitted trapezoidal rule", Includes filtering of Ill-posed 
Initial conditions, automatic stepslze selection, and automatic selection of 
Newton-Jacobl or Newton Iteration for convergence to achieve maximum computa- 
tional efficiency while observing a prescribed error tolerance. The new algo- 
rithm was found to compare favorably with LSODE on two representative test 
problems drawn from combustion kinetics. 


INTRODUCTION 

The problem of economical Integration of the coupled, nonlinear ordinary 
differential equations describing exothermic, homogeneous gas phase combustion 
reaction kinetics has not yet been optimally solved. Multlstep methods, based 
on traditional explicit numerical Integration schemes such as the fourth-order 
Runge-Kutta method, dominated the thinking of early workers performing single- 
point calculations modeling one-dimensional combustion processes In shock tubes 
and rocket nozzles. In the 1960's, considerable progress was made with respect 
to reliability and efficiency of single-point calculations, due to the Intro- 
duction of Implicit methods by Tyson (ref. 1), and Treanor's Introduction of 
locally exact solutions (exponential functions) to extend the convergence ra- 
dius of explicit Runge-Kutta methods (ref. 2). Lomax and Bailey (ref. 3) ex- 
plored the combined use of explicit and Implicit methods based on low-order 
polynomial approximations as an alternative to Treanor's method. 

In the 1 970 1 s , Treanor's algorithm was utilized by Dimitrov to model hy- 
drogen-air combustion (ref. 4), and Tyson's and Lomax and Bailey's Ideas were 
Incorporated Into production codes (refs. 5,6). The first Implementation 
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of Gear's variable-order, variable-step, backward-difference formulas 
(refs. 7,8) was produced by Lawrence Livermore Laboratory (ref. 9). 

At the same time, Increasing attention was focused on multipoint calcula- 
tions, required for modeling of multidimensional reactive flows (refs. 10-12). 
This Increasing Interest In modeling of furnaces and combustors resulted In a 
major reassessment of available algorithms. It was recognized that all pro- 
posed techniques for modeling of complex chemically reacting f low--whether 
Eulerian space discretization with operator splitting for treating coupled 
phenomena, or Lagranglan mass or vortlclty discretization with coupled phenom- 
ena treated by the method of fractional steps--all require a very stable but 
only moderately accurate homogeneous batch chemistry Integrator (refs. 10-12). 
The same requirements exist for the single-point, stochastic simulation of 
turbulent mixing-influenced. Inhomogeneous gas-phase continuous combustion 
systems (refs. 13-16). 

At the present time, the latest of a series of Lawrence Livermore 
Laboratory's Implementations of the Gear algorithms, LSODE (refs. 17,18) is 
considered to be the best available, state-of-the-art code for solving stiff 
systems of ODE's. However, It Is recognized by combustion device modelers 
that LSODE Is not fast enough for economical calculations of multidimensional 
reacting flowflelds. 

While LSODE Is regarded as the best available "packaged" code for solving 
an arbitrary system of ODE's, It may be possible to construct a superior method 
for solving a uniquely specified system of ODE's. A clue to a better method 
for the present problem than the variable order-variable stepslze algorithm 
utilized by LSODE Is given by Lambert (ref. 9, p. 18): 

"For stiff equations, for which a typical solution Is a rapidly decaying 
exponential, the error In Interpolating such functions by a polynomial of 
high degree Is very large, and correspondingly we find that linear multi- 
step methods of high stepnumber cannot be used successfully for such 
problems . " 

The present paper describes the development of an algorithm, denoted 
CREK1D (for Combustion Reaction Klnetlcs-One-Dlmenslonal ) , which has been spe- 
cifically designed to Integrate chemical kinetic rate equations, and which 
compares favorably with LSODE. 


GOVERNING ALGEBRAIC AND DIFFERENTIAL EQUATIONS 


The equations for adiabatic, homogenous gas-phase chemical reaction ki- 
netics at constant pressure are given by: 


= VVT). ^ = l.N (1) 
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where 


(3a) 


(3b) 

In equations (1) to (3), Is the mole number of the ith species 
(1 = 1,N); T Is the temperature; p Is the mass density; o^j and 
are the stoichiometric coefficients of species 1 (1 = 1,N) In reaction 0 
j(j = 1,J) as a reactant and as a product species, respectively; Aj, A_j, Bj, 
B_j, Ej, and E_j are constants In the modified Arrhenius rate expres- J 
slons for Rj and R_j, which In turn are the forward and reverse rates 
of the jth reaction (j = 1,3) as a reactant and as a product species, respec- 
tively. N and J are, respectively, the total number of distinct chemical 
species In the gas mixture and the total number of distinct elementary 
reactions In the reaction mechanism. 

For adiabatic constant-pressure chemical reaction, the following enthalpy 
conservation equation constitutes an algebraic constraint on equations (1) to 
(3): 


°1 " a k1 

Rj = AjT j exp(-Ej/RT) IT (p« k ) KJ 

k — 1 
and 

B , N a' ' 

R _j = A _j T exp(-E__j/RT) TT (po k ) KJ 


N _ 

Xh.d, = H = const. (4) 

1=1 1 1 0 

where Ti^ Is the molal-speclf 1c enthalpy of species 1, and H 0 Is the 
mass-specific enthalpy of the mixture. 

The mass density p In equations (2) and (3) Is determined by the equa- 
tion of state for an Ideal gas. 


P = P/(RTo m ) (5) 

where P Is the absolute pressure, R the universal gas constant, and 

°m = °1 (6) 

Is the reciprocal mean molar mass of the mixture. 


PHYSICAL AND COMPUTATIONAL SCENARIOS 

A typical constant pressure batch combustion problem consists of three 
distinctly different chemical-physical regimes: Induction, heat release and 

equilibration (see fig. 1). To accommodate the widely different characteris- 
tics of these regimes, a two-part algorithm was developed as follows. 
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During Induction and early heat-release, the species equations are domi- 
nated by positive time constants, and the temperature also exhibits a positive 
time constant. Since very small steps are required for Integrating unstable 
equations, a simple predictor-corrector scheme with functional Iteration as- 
sures the least computational work possible. However, during late heat release 
and equilibration, when the temperature and species equations exhibit negative 
time constants, large stepslzes can be used, so Newton-Raphson Iteration with 
calculation of the full Jacobian matrix Is the optimal convergence method. 

Although the governing ordinary differential equations (ODE's) are stable 
during late heat release and equilibration, they are characterized by widely 
differing time constants. This behavior— termed "stiff" by Curtiss and 
Hlrschfelder (ref. 20)--causes classical Integration techniques (such as the 
popular explicit Runge-Kutta method) to use very small steplengths, thereby 
resulting In excessive computational work (refs. 19, 21-23). 

The equilibration process does not have a clearly defined termination, 
due to the asymptotic nature of the approach to the chemical equilibrium state. 
Since equilibrium values of temperature and species concentration can be de- 
termined a priori by a Gibbs function minimization scheme (refs. 11,24), the 
end of the equilibration period can be defined as the time at which all of the 
mole numbers and the temperature are within (say) one percent of their chemi- 
cal equilibrium values. 


CANDIDATE ALGORITHMS 

A number of single-step algorithms were considered. Including: 

1. Implicit Euler rule (refs. 1, 3, 19, 25, and 26) 

2. Implicit midpoint rule (refs. 27, and 28) 

3. Trapezoidal rule (ref. 19) 

4. Exponential-fitted trapezoidal rule: "Llnlger-Wllloughby No. 1" 

(refs. 12, 25, 26, and 29-31) 

5. Trapezoidal rule with end correction: cubic spline, or "Obreschkoff 

(2,2)" (refs. 25, 28, and 29) 

6. Exponential-fitted trapezoidal rule with end correction: "Llnlger- 

Wllloughby No. 3" (refs. 25, 26, and 29) 

7. Local ly-exact or exponential solutions (refs. 19, 25, and 26). 

Variations and combinations of the seven basic algorithms were explored; 
for example. Young and Boris' SAIM (Selected Asymptotic Integration Method), 
which utilizes a combination of local ly-exact ("asymptotic") solutions and an 
explicit Euler-trapezoidal rule predictor-corrector method (refs. 10, 32, 33). 

Of the seven algorithms and their variants tested. It was found that the 
most useful algorithms were variations of "Llnlger-Wllloughby No. 1," the 
exponential-fitted trapezoidal rule. 
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METHOD FORMULATION: THE EXPONENTIAL-FITTED TRAPEZOIDAL RULE 

Consider an approximate solution to equation (1) based on a variation of 
the trapezoidal rule, the "tuneable trapezoid" (refs. 12, 25, and 26) or 
"theta-method" (ref. 28), as follows: 

°1 ,ntl ■ * "(Vt.nvi *«’ - 1 - 1 -* (7) 

where, o\ n Is the approximation to the exact solution to equation (1) at 
the current time, t n , h Is the time steplength (= t n+ -|- t n ), f^ >n = fl(<*k,n* T n)* 
and T n , the temperature at t n . Is obtained from equation (4) — see 
appendix A for details. Also, Is a degree-of-lmpllcltness or "tuning" 
factor such that = 0 recovers the explicit Euler approximation, U^ = 1 
recovers the Implicit Euler approximation, and U^ = 1/2 recovers the trape- 
zoidal rule (or modified Euler method). 

Following the Ideas of Llnlger and Willoughby (ref. 29) and of Brandon 
(refs. 30 and 31), we Introduce the concept of "exponential-fitting" the para- 
meter to a local ly-exact solution of equation (1). Assume a locally- 

linearized form of the rate equations, equation (1): 

f 1 = (f 1,n " Vl.n* * Vi (8a) 

or 

cTT = (f 1,n " Vl,n* + Vi (8b) 

where the choice of , a suitable linearization constant, Is discussed 
In detail In the next section. Integration of equation (8b) gives the follow- 
ing result 


o. , = o. + hf, 
1,n+l 1,n 1,n 


exp(e^h)-l’ 

e^h 


1 = l,N 


which Mlranker (ref. 26) calls the "filtered Euler" approximation. 

Application of equation (8a) to the entire step of length h gives 


(9) 


f. , = (f. - 0.0. ) + 9.0. . 

1,n+l v 1,n 1 1,n' 1 1,n+l 


(10) 


Substituting equation (10) for the term f^ .j In equation (7), eliminating 
o^ between equations (7) and (9), and solving explicitly for 


1 1 ■ i 

U 1 ~ e^h + l-exp(e^h) 


which relates the tuning factor In equation (7) to the linearization 
constant In equation (8). Equation (11) Is graphed In figure 2. 
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However, In order to maintain the absolute A-stabllity of equation (7) 
(i.e., cm n+i remains bounded as h Is Increased Indefinitely), It Is nec- 
essary to’restrlct to the Interval (0.5, 1.0) (refs. 19, 25, 26, and 
28-31). Thus, equation (7) Is made to default to the second-order-accurate 
trapezoidal rule whenever Is greater than zero. However, whenever 
ei Is negative, equations (7) and (11) together are equivalent to the 
locally-exact or exponential solution, equation (9). Brandon (refs. 30 and 
31) has shown that the equivalent polynomial accuracy of equation (9) Is typi- 
cally of order six to eight. 

Thus, equations (7) and (11), with the restriction (0.5 < < 1), con- 

stitute an exponential-fitted trapezoidal rule, a method which Is A-stable, 
and has a polynomial-order accuracy of at least two, and as great as slx-to- 
elght. It Is Interesting to note that, using Llnlger and Willoughby's termi- 
nology, Treanor's method (ref. 2) would be called an "exponential-fitted, 
fourth-order Runge-Kutta" method. 


LINEARIZATION OF THE RATE EQUATIONS 

There are at least two distinct ways of Interest to determine the line- 
arization constants In equations (8) to (11). First, the traditional 

"chemist's" approach, which may be termed a formal linearization. Is reviewed. 


Formal Linearization of the Rate Equations: "L-Formulatlon" 

Equation (2) may be expressed as a difference between two positive- 
definite terms as follows (refs. 10, 32, 33) 


fi » Ql - Dl 

(12) 

where 


Q 1 E p 1 (a 1j R -j * “ijV 

(13) 

D 1 = P 1 + a 1j R -j } 

(14) 

The terms and represent the gross rates of production and 

consumption of specles-i, respectively, due to the contributions of all J 
forward and reverse reactions. The objective of this decomposition Is to en- 
able factorization of the mole number from the destruction term: 

V L i*i 

(15) 

where L^ , obtained simply by dividing by a \ , Is 

given by 



With the notation of equations (12) to (16), equation (1) may be rewritten 


as: 




Vi 


(17) 


Equation (17) Is now formally linear In 01 , and It can be seen by Inspection 
that (-Li) Is the appropriate linearization constant ©i . With replaced 
by (-Li), and fi replaced by (Qi - Lioi), equation (8a) may be rewritten 
In the form 


°1 ,n+l = [Q 1,n /L 1,n ] + [o 1,n 


Q 1,n /L 1,n ] exp( " L 1,n h) 


(18) 


Since the loss coefficients (Li, 1 = 1,N) are all positive-definite, 
the first terms In brackets in equation (18) may be regarded as the local as- 
ymptotic solution (refs. 10, 32), which represents the large-tlmestep value of 
the locally exact solution, equation (18). 


Inspection of equation (18) reveals that the coefficients (Li , 1 = 1,N) 
represent Inverse characteristic time constants for the L-formulated rate 
equations; that Is, an accurate approximate solution according to equation (18) 
would have to be resolved on a steplength of order l/(L^) max . 


Functional Linearization of the Rate Equations: "Z-Formulatlon" 

Integration of equation (17) to yield equation (18) required that Qi 
and Li In equation (17) be constant over the tlmestep h. It very often 
happens that this condition does not exist. Therefore, a more careful exami- 
nation of the actual coupling between variables In equation (1) requires a 
functional linearization of the rate equations. 

Equation (10) can be solved explicitly for to give 


e 


1 


f, , - f< 

LlELLI l.n E z 

°1 ,n+l tf 1 ,n 


(19) 


The parameter 1\ Is termed by Brandon the "diagonal transition matrix" 
or the "state variable differential" (ref. 30). Also, may be recognized 

as the average value over the interval h of the ratio of second-to first-time 
derivatives of the mole number 

d 2 o,/dt 2 df ,/dt 

“do^/dF = “Tf (20 > 


Therefore, Zi Is a measure of the local curvature of the trace of 
with time, as for example In figure 1. The parameter Zi Is also 
equivalent to Treanor's parameter Pn (ref. 2). 

An Important relationship which follows from the locally exact solution 
Is obtained by substituting equation (9) Into equation (10) (with Q\ = Zi n ) 
to give ’ 
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(21) 


fl,n+l = U,n exp (Z i>n h) 

Equation 21 emphasizes the role of the state variable differential Z^ as a 
"diagonal transition matrix" (refs. 30, 31). Thus may be thought of as 
a finite-difference analog of the eigenvalues of the Jacobian matrix of 
equations (1) and (2). 


Comparison of Land Z-Formulatlons 

Because the loss coefficients L^ , equation (16), are all positive- 
definite, the Integrated rate equations, equations (7) to (11), are always 
stable when = -L^. However, they are not always accurate, because 
strong coupling between variables often results In observed positive time con- 
stants Z] for traces of o\ versus t1me--see figure 1. 

In contrast to the formally-linearized equations (e^ = -L^), the 
functionally-linearized equations (e^ = Z^) are always accurate, but are 
stable only when Z^ Is negative; when any species exhibits a positive 
l\, both the physical equation and the corresponding approximate solution, 
equations 7 to 11 with = Z \ , are unstable. As a consequence, when 
the rate equations are dominated by positive Z^'s, as occurs during Induc- 
tion and early heat release. It Is necessary to take small steplengths of the 
order 1 / I Z^ | max • 

The use of the Z-formulated equations may result In extremely small step- 

lengths (of size 1/Uilmax) be1ng taken 1f the in1t1a1 va1ues for 
o^(1 = 1 ,N) are Ill-posed. These non-physical Initial conditions may arise, 
for example In multidimensional modeling because of the averaging of mole num- 
bers over adjacent grid nodes. In this case, use of 0^ = -L^ with a 
time steplength of order l/|L^| max Is desirable. On the first call 
to CREK1D, the L-formulated equations are solved over one tlmestep of length 
1 / | L 1 | max to t1 lter the Initial cond1t1ons--that Is, to remove spuri- 
ous transients of time scale 1 / 1 | ma x and to provide physically 
meaningful Initial values. 

Another difficulty with the use of the Z-formulated equations arises 
whenever any species Is In "quasi-steady state"; that Is, when the time trace 
of a\ closely approaches an asymptote or passes through an extremum. In 
either case, 1\ becomes singular due to division by zero or a small number 
(see eq. 19). When this condition occurs — that Is, when f^ ~ 0 and/or 
(<m n+1 - oi >n ) ~ 0--Qi and (eq. 17) are locally only slowly 

varying wlth’tlme, so that use of the L-formulated equations (l.e., 0^ = 

-L^), which are now both stable and accurate, Is appropriate. 


SOLVING THE INTEGRATED RATE EQUATIONS 

A fundamental question which must be addressed Is whether or not to use 
Newton-Raphson Iteration to converge the Implicit equations, equations (7) to 
(11), together with the enthalpy conservation constraint, equation (4). The 
tradeoffs which must be considered Involve accuracy, convergence radius and 
rate, and computational work. 


8 



Newton-Raphson (NR) Iteration Is attractive because It converges quadrat- 
Ically and has an Infinite convergence radius. However, NR Iteration requires 
frequent evaluation of the Jacobian matrix and Its Inverse, usually done by 
either Gaussian elimination or (equivalently) LU-decomposItlon and back- 
substitution (ref. 34). 

On the other hand, functional Iteration techniques such as nonlinear 
Gauss-Seldel , Jacobi, or Jacobl-Newton (JN) Iteration (ref. 35) do not require 
evaluation of the Jacobian matrix or Its Inverse, but have severely restricted 
convergence radii. In addition, for these methods, the convergence rate Is 
only linear, or at best, super 11near--better than linear but not quite 
quadrat1c--1n the case of JN Iteration. 

For reasons presented In the section Physical and Computational Scenarios, 
JN Iteration Is used during Induction and early heat release; however, during 
late heat release and equilibration, NR Iteration Is used. Details of both JN 
and NR Iteration methods are given In appendix A. 


Approximations for State Variable Differentials 

The "tuning factors" In equation (7) are given by equation (11) 

with representing the corresponding state variable differentials 
Z'j f p . To minimize computational work, the Z^ >n 's are evaluated only 
once per step--at the beginning of the time step, using equation (19). How- 
ever, since o^n+1 an d ^l.n+l are n °t known at the start of the step, 
an approximation has to be developed for Z^ >n . This Is done simply by using 
values from the previous step, so 


Z 



f 1 ,n-l 

-a, , 

1 ,n-l 


Selecting the Iteration Technique 


( 22 ) 


CREK1D automatically selects the Iteration scheme (JN or NR) to be used 
for solving equation (7). During Induction and heat release, when small step- 
lengths are required for solution stability (refs. 21, 23), JN Iteration Is 
used to minimize computational work. During late heat release and equilibra- 
tion when the differential equations are more stable and larger steplengths 
can be used (refs. 21, 23), NR Iteration Is preferred since It has a much 
larger radius of convergence than JN Iteration. The regime Identification 
test exploits the fact that during extrema or equilibration ("quasi-steady 
state": QSS) many reactions achieve a condition In which the forward and re- 

verse reactions are large but with vanishingly small differences (refs. 36, 

37). The actual test employed at the beginning of the time step Is 

1^1 < 10" 3 (Q i + D 1 ) (23) 

where, and D^ are the production and destruction rates, respectively, 
of species 1 (see eqs. 12 to 14). If any two are In QSS— that Is, satisfies 
equation (23)— NR Iteration Is used for the step. If fewer than two species 
satisfy equation (23) JN Iteration Is used for the step. 
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Non-Physical Initial Conditions 

As discussed In the section Comparison of Land Z-formulatlons, nonphysical 
Initial conditions may result In extremely small steplengths being taken. To 
filter the Initial conditions — that Is, to provide physically meaningful Ini- 
tial mole numbers and net species production rates--the L-formulated equations 
are solved over one tlmestep. On the first call to CREK1D, It uses this for- 
mulation over a tlmestep of length h-| given by 

h, - (20 


'1 


max 

1 


-1 


The predictor-corrector algorithm uses equation (18) as the predictor 


o\°\ = ^(0) * h-j f ( 0) 


"1 - exp(-L 1 (0)h 1 )“ 


(25) 


An Implicit Euler corrector equation Is then Iterated to convergence 

(m+1) , n . . f (m+l) 


(26) 


In these equations oi(0) are the Initial values, f 1 ( 0) = fi(ok(0), 

T(0)), T(0) Is the Initial temperature, and the subscript 1 Is used to Indi- 
cate that this Is the first step. Using equations (17) and (A12) 

together with the approximations Q^ m j^= j and l!| j ^ ], equation (26) 

can be rewritten to provide the following expression for the log-variable 

, (m) 

corrections Alog -j 

(m). i_ £ (m) 


tn\ V'"/. h f''"' 

Alna (m) °1 (0) I q 1.1* h l f 1,1 . ,. 1N 

Alog °1,1 = ( m ) + h D (m) ’ 

° 1,1 + Vl.l 


(27) 


Equation (27) Is Iterated until converged (see appendix A). If convergence Is 
not achieved after 10 Iterations, the steplength Is halved and the step re- 
tried. If convergence Is obtained after M (M < 10) Iterations, the step Is 
accepted as successful, and the solution for Is updated 


* 1,1 


,(M). 

1.V 


= <*;;; 1 = l,N 


(28) 


and the temperature T-j Is obtained by a single Newton-Raphson Iteration 


H 


T 1 = T(0) + 


N 

E * 

1=1 




E ° i , iV r<0)) 

1=1 


( 29 ) 


where c^ Is the constant-pressure molal specific-heat of species 1. 
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Species In Quasi-Steady State 

Whenever the reaction rate for any species 1 satisfies equation (23), 
that species Is considered to be In "quasi-steady state" (QSS). For reasons 
presented In the section Solving the Integrated Rate Equations, the 
L-formulated equations (e^ = -L^) are employed for all species In quasi- 
steady state. 


ACCURACY, STABILITY, CONVERGENCE AND STEPLENGTH CONTROL 

It Is particularly useful to define two timescales needed for steplength 
control: "h acc y," the estimated maximum stepslze to stay within a pre- 

scribed local truncation error tolerance e, and "hiter»" the maximum step- 
size which will admit efficient convergence (refs. 38, 39). 


Local Truncation Error (lte) 


In the nonstiff regime where JN Iteration Is used, many species equations 
have positive time constants. Such species, which we term unstable, have 
values of > 0. In the present version of the code, whenever 1\ > 0, 
the "tuning factor" (eq. 7) Is set equal to 0.5, so that with the fol- 
lowing relations deduced from equation (11) (with = Z^), 


exp(Z^h) - 1 1 

Z^h = 1 - U^h 


1 


1 - 


\ z i h 


(30) 


the predictor equation (A21) for unstable species becomes 


a 


( 0 ) 

1 ,n+l 


a. + 

1,n 





(31) 


where Z^ >n Is determined by equation (22). The corrector equation for 
unstable species Is equation (7) with = 0.5--that Is, the trapezoidal 
rule 


a 


(m+1 ) 
1 ,n+l 



(32) 


It can be shown that equation (31) has approximately the same lte as a 
second-order Adams-Bashforth predictor, so the lte for unstable species can be 
conservatively estimated as 


(lte) 1 


1 

6 


°1 ,n+l 


( 0 ) 

1 , 


(33) 


where Is the converged solution. 

Species equations that are stable are not limited to = 0.5, so 
equation (33) does not apply for such species. In fact, no such estimate as 
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equation (33) can be made for the present method because It Is not a poly- 
nomial method (refs. 30 and 31). Nevertheless,, equation (33) Is used even 
though It seriously overestimates the lte for stable species. 

For all species, stable and unstable, an average weighted lte estimate Is 
made by 


lte 



(34) 


and the projected next steplength which would satisfy the user-specified local 
error tolerance e Is calculated from 


h = h(e/lte) 1/3 (35) 

accy 

In the computationally stiff regime where more of the species are stable, 
and where use of the NR Iteration permits much larger steplengths, equations 
(33) to (35) are far too conservative. In fact, In this regime, the local 
truncation error for stable species may decrease with an Increase In the step- 
length. Therefore, a different estimate of the lte Is used In the stiff 
regime. The estimate Is based on the norm of the difference between the con- 
verged solution and that obtained after the first Newton Iteration. This 
estimate Is actually a measure of the local linearization error due to the use 
of NR Iteration, and not that of the local truncation error. The relations 
used In the stiff regime In place of equations (34) and (35) are 


and 


lte = 



(1) 


i N 

i v 

°1.n+l ^l.n+l 

1 2 

n 1=1 | 

max(cr. , a. , ) 

v 1 ,n’ 1 ,n+r 



h aC cy - h(«/lte) 


1/3 


(36) 


(37) 


In equations (36) and (37), the factor (1/3) must be regarded as strictly 
empirical . 


Convergence Control 

During JN and NR Iteration the rate of convergence Is monitored, both to 
detect divergence and to optimize the rate of Iterative convergence. 

Following Shamplne (refs. 38, 39) and Pratt (ref. 40), the convergence rate 
R CO nv defined by 
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R 


conv 


Vs ^ 


(38) 


where (m-1) and (m) denote the Iteration numbers. Note that at least two 
Iterations are required to define R conv ; If convergence Is obtained after 
the first Iteration, a default value of R C onv = 0-1 for Iteration and 
0.05 for NR Iteration Is used. As with the convergence test, the summations 
In equation 38 Include only species with mole numbers greater than lO -2 ^. 


If Rconv greater than one, the solution Is diverging and the step 
must be attempted with a decreased steplength. In CREK1D, a more conservative 
test R CO nv > 0*8 Is used to detect divergence. If convergence Is not ob- 
tained or If divergence Is detected, the steplength Is decreased as discussed 
In the next section. 


Steplength Control 

During both JN and NR Iteration the convergence rate, R conv , used to 
control the steplength. If corrector convergence Is not obtained after UMAX 
Iterations, where UMAX Is the user-supplied maximum number of corrector 
Iterations to be attempted, or If divergence Is detected (1.e. t R conv > 0.8), 
the steplength Is decreased. The new steplength, h 1 , Is calculated as follows 

h' = h min jo. 5, max(0.1, 0 . 5 / R con v )^ (39) 

and the step retried with the decreased steplength. 


After corrector convergence, the steplength Is adjusted up or down as 
necessary to keep R CO nv the range (0.4, 0.5). An estimate, h^ er , Is 
made of the steplength that would result In the desired convergence rate: 


h(0.4/R ) 

' conv 7 


’Iter 


= h 


h(0 - 5/R conv> 


1/2 


1/2 


R < 0.4 
conv 

0.4 < R < 0.5 
conv — 

R rnn v/ > 0 - 5 

conv 


(40) 


In addition, the maximum permissible steplength, h accy (given by eq. 35 or 
37) that would result In a local truncation error equal to the user-specified 
value for the local relative error, e, allowed per step, is calculated. 


The steplength, h 1 , to be attempted for the next step is then taken to be 
the minimum of h^er an( * ^accyJ however, the ratio h'/h Is restricted 
to be no larger than 10. Hence, h' Is given by 

h' = min (h^ er , h accy , lOh) (41) 
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To Inhibit oscillations of the steplength, h (due to repeated failures of 
the convergence test after Increases In h), an Increase In h Is allowed only 
If convergence was obtained In the previous step without a reduction In h. 
Following any step for which corrector convergence could not be obtained or 
divergence was detected (and so h was reduced), the new steplength h' Is 
calculated from 


h' = min (h, h' given by eq. 41) (42) 


COMPUTATIONAL STRATEGY 

The objective of the calculation Is to efficiently and reliably solve the 
following Initial value problem: given a set of Initial conditions (o^; 1 = 

l.N), T and P (constant), find the values of (1 = 1,N) and T at the 
end of a prescribed time Interval At. The overall steps taken to accomplish 
this are: 

a. Filter the Initial conditions: evaluate the f^'s; choose a time- 

step h-| = l/(L^) max ; use the "filtered Euler" predictor, equation (9) 

with = -Li (1 = l.N), and converge equations (A19) and (A20) by JN 
Iteration, with all Ui's set to unity (Implicit Euler corrector). 

b. Adjust the tlmestep according to the observed R conv from step a. 
Use an explicit Euler predictor (eq. 9) with =0(1 = l.N); use JN 
Iteration, equations (A19) and (A20), with all Ui's set equal to one-half 
(trapezoidal rule). 

c. Calculate R conv , ^accy an d ^1 ter I choose JN or NR Iteration, 
according to whether or not two or more species are In QSS. Adjust the time- 
step as described In the section Steplength Control. 

d. Evaluate Zi's from equation (22), Ui's from equation (11) and 
Iterate equation (17) until converged. 

e. Repeat step d to the end of the prescribed Interval At; return 
solution. 

The strategy outlined above. In conjunction with the robust Integration algo- 
rithm, leads to minimum computational time with acceptable accuracy. 


COMPUTATIONAL TACTICS 

The computational work of evaluating logarithms and exponentials may be 
avoided by judicious use of approximating functions. For example, the term 
(e x - l)/x In equation (9) Is evaluated In the code by means of a (2,2) di- 
agonal Pade' (rational function) approximation for e x , 


e 


x 

( 2 , 2 ) 



(43) 
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which results In 


r 



] 

1 - x(l/2 - x/1 2) ' 


x < 0 


(44) 


Similarly, the "tuning factor" of equation (11) can be evaluated with 
sufficient accuracy by means of an approximation which has no singularity 


1 

x 


+ 




(45) 


The right-hand side of equation (44) requires only five operations to 
evaluate, and does not exhibit the singularity at x = 0 of the exact left- 
hand side. Equation (44) Is used only for the predictor equation (9). 

Equation (45) Is graphed In figure 2 ; note that equation (45) defaults to 
= 0.5 (the trapezoidal rule) for x > 0. 

Another significant reduction In computational work was achieved by lo- 
cally linearizing the fifth-order polynomial approximations used In evaluating 
the species-1 enthalpy h^ and specific heat capacities c^, which appear 
In equations (All) and (A20). 

Finally, It should be noted that, while log-variable corrections are used 
In the algorithm, (appendix A), evaluation of logarithms of the variables 
themselves Is avoided by use of the approximate equations (A12). The only 
evaluation of exponentials In the code are those necessarily required for the 
Arrhenius rate expressions In equation (3) and those required for the tuning 
factor, equation (45). 


PERFORMANCE OF CREK1 D— COMPARISON WITH LSODE 

A preliminary version of CREK10 (ref. 40) has been tested by Radhakrlshnan 
(refs. 21 to 23) against LSODE (refs. 17,18) on two test problems drawn from 
combustion kinetics. Both problems described adiabatic, constant pressure, 
transient batch chemical reaction and Included all three regimes of 
combust1on--1nduct1on, heat release and equilibration. 

Test problem 1, Illustrated In figure 1, described the Ignition and sub- 
sequent combustion of a mixture of 33 percent carbon monoxide and 67 percent 
hydrogen with 100 percent theoretical air, at a pressure of ten atmospheres 
and 1000 K Initial temperature. It consisted of 12 reactions Involving 11 
species. Test problem 2 described the Ignition and subsequent combustion of a 
stoichiometric mixture of hydrogen and air at a pressure of two atmospheres 
and 1500 K Initial temperature. It Involved 30 reactions among 15 species. 

Both test problems were solved over a time period of 1 ms In order to obtain 
near-equilibration of all species and the temperature. 

In applying LSODE to the problem of solving chemical kinetic rate equa- 
tions, two different methods (A and B) for calculating the temperature were 
attempted. In method A (LS0DE-A), the temperature was calculated from the 
mole numbers and the Initial mixture enthalpy using the algebraic enthalpy 
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conservation equation (4) and a Newton-Raphson Iteration technique. In method 
B (LSODE-B) the temperature was evaluated by Integrating Its time-derivative, 
obtained by differentiating equation 4 with respect to time 


dT 

dt = 


N 

- 13 h^f^ 
1=1 1 1 



(46) 


In applying the present version of CREK1D to the two test problems dis- 
cussed above, we have adopted the procedure described by Radhakrlshnan 
(ref. 23) and summarized here. A typical computational run consisted of Ini- 
tializing the species mole numbers, temperature and CPU time. The Integrator 
was then called with values for the necessary Input parameters. Including the 
local error tolerance, e, required of the numerical solution and the elapsed 
time (= 1 ms for both problems) at which the Integration was to be terminated. 
On return from the Integrator, the total computer time (CPU) required to solve 
the problem was calculated. 


Figures 3 and 4 present the computational work (expressed as the CPU time 
In seconds required on the NASA Lewis Research Center's IBM 370/3033 computer) 
required by the single-precision codes CREK1D and LSODE, plotted against the 
local relative error tolerance, e. 


Figures 3 and 4 show that CREK1D compares favorably with LSODE for large 
values of the relative error tolerance, e. But for small values of e 
CREK1D Is slower than LSODE. However, use of low values of e Is wasteful 
because of uncertainties In reaction rate coefficients (ref. 10). In addition, 
the proposed purpose of CREK1D Is to perform multipoint calculations of chemi- 
cally reacting flows by coupling It with a hydrodynamic equation solver. These 
solvers are at best accurate to within a few percent, so generation of highly 
accurate chemical kinetic solutions Is wasteful (ref. 10). 

The solution of the coupled hydrodynamic-reaction rate equations requires 
the solution of the reaction rate equations at several thousand grid points 
for relatively short periods of time (ref. 10). Hydrodynamic codes also have 
large storage requirements. Hence, reaction rate Integrators with both a small 
storage requirement and a low Initialization (start-up) time are needed. The 
storage and start-up time requirements of CREK1D have been shown to be signif- 
icantly less than those required by LSODE (ref. 24). These factors make CREK1D 
more attractive than LSODE for multipoint calculations. 

To further explore the differences In computational work required by 
CREK1D and LSODE, we present In figures 5 and 6 plots of the steplength suc- 
cessfully used by these codes through the course of each problem. For test 
problem 1 (fig. 5) and for test problem 2 at long times (fig. 6), the step- 
lengths selected by CREK1D are comparable to those selected by LSODE. How- 
ever, at early times (t ;$ 5 ys) for problem 2, CREK1D uses much smaller 
steplengths than LSODE and Is hence slower. 
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The present version of CREK1D will be available for testing from COSMIC, 
University of Georgia, Athens, Georgia 30602. Details of the subprograms 
included In CREK1D are given In appendix B. In appendix C, we present a 
sample Input data set and the resulting output. These were generated on the 
NASA Lewis Research Center's IBM 370/3033 using single-precision accuracy. 


CONCLUSIONS 

A major conclusion of the present work Is that the L-formulated equations, 
widely used at present for modeling both atmospheric chemistry (refs. 32 and 
33) and combustion processes (ref. 10), should be employed In the Initial In- 
duction period, when some mole numbers may be very small, and whenever any 
species Is In quasi steady-state. At all other times, use of the Z-formulated 
equations completely obviates the need to use asymptotic or quasi-steady-state 
assumptions to resolve the near-equilibrium stiffness problem uniquely associ- 
ated with the L-formulated equations (refs. 10, 27, 32, and 33). For atmos- 
pheric chemistry problems In particular, the present algorithm without the 
enthalpy constraint should constitute a very fast and robust method for calcu- 
lation of Isothermal, batch homogeneous gas-phase kinetics. 

The success of the present algorithm stems from recognition of the fact 
that the approximating equations describing the Induction processes are 
stability-limited, whereas the corresponding post-induction equations are 
accuracy-limited. In the former case, low-accuracy methods with the least 
possible computational work per Iteration are Indicated due to the necessarily 
restricted stepslze, as pointed out by Young and Boris (ref. 10). However, In 
the post-induction processes, where Inherent stability Is not a problem, step- 
sizes may be as large as the accuracy and convergence radius of the approxi- 
mating equations permit. The use of exponentials as approximating functions 
satisfies the accuracy requirement, and Newton-Raphson Iteration gives an 
infinite convergence radius. The added computational work of Jacobian evalua- 
tion and matrix Inversion required for Newton-Raphson Iteration Is offset by 
very large stepslzes and quadratic convergence (refs. 19, 25, and 26). 


NOMENCLATURE 

Aj » A-j preexponentlal constants In forward and reverse rate equations 

for reaction j (eq. 3), units depend on reaction type 

B j , B_j temperature exponent In forward and reverse rate constants for 
reaction j (eq. 3) 

CPU total CPU time required on IBM 370/3033, s 

CpI constant-pressure molal-speclf 1c heat of species 1, J/kmole K 

Dj rate of destruction of species 1 (eq. 12), kmole 1/kg mixture s 

E j , E_j activation energy In forward and reverse rate equations for reac- 
tion j (eq. 3), cal/mole 

fl net rate of formation of species 1 (eq. 2), kmole 1/kg mixture s 

H 0 Initial mixture mass-specific enthalpy, J/kg 
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haccy 

hlter 


/~s_/ 

hi 

UMAX 

J 


Li 


1 te 1 

N 

P 

Ql 

R 

^conv 
R j* R -j 

T 

t 

Ul 

*1 

I I I 

<*1j» a 1j 


®1 

p 

°1 

°m 


steplength used by Integrator, s 

i 

estimated steplength that would result In solution with desired 
accuracy (eqs. 35 and 37), s 

estimated steplength that would result In a convergence rate In 
the range (0.4, 0.5) (eq. 40), s 

molal-speclf 1c enthalpy of species 1, 3/kmole 

maximum number of corrector Iterations to be attempted by the 
Integrator 

total number of distinct elementary reactions In reaction 
mechanism 

loss coefficient for species 1, Inverse of characteristic time 
constant for species 1 (eq. 15), 1/s 

local truncation error for species 1 

total number of distinct chemical species In the gas mixture 
absolute pressure, N/m 3 

production rate of species 1 (eq. 12), kmole 1/kg mixture s 
universal gas constant, 8314.3 J/mole K (1.9872 cal/mole K) 
Iteration convergence rate (eq. 38) 

molar forward and reverse rates per unit volume for reaction j 
(eq. 3), kmole/m 3 s 

temperature, K 

time, s 

degree of Implicitness or tuning factor for species 1 (eq. 7) 

state variable differential for species 1 (eq. 19) 

stoichiometric coefficients for species 1 In forward and 
reverse reaction j (eq- 3), number of kmoles 1 In elementary 
reaction j as a reactant and as a product, respectively 

local relative error tolerance 

suitable linearization constant for species 1 (eq. 8) 
mixture mass density, kg/m 3 
mole number of species 1, kmole 1/kg mixture 
reciprocal of mixture mean molar mass (eq. 6), kmole/kg 
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APPENDIX A 


OUTLINE OF NR AND JN ITERATION TECHNIQUES 
Newton-Raphson Iteration 

A Newton-Raphson functional Is defined from equation (7) by 

(m) 


F (m) _ ”1.n+1 °1.n 
1,n+l ~ ^h 


1-U, 


- f 


(m) 


1,n 1,n+l 


(Al) 


where m Is the Iteration number, fj m * +1 = IV^n+l’ T n?l^’ and 

= 1»N) and T n+1 are a PP rox1mat1ons to the exact values 

<> k (t n+1 ) and T(t n+1 ), respectively. 

The Newton-Raphson functional for temperature Is determined from equation (4), 


c( m ) _ y rt (m) r / T (mh H 
F T,n+l - "l.n+l h 1 (T n+l } - H 0 

where hi (Tn?j) Is the molal-speclf 1c enthalpy of species 1 at 

temperature Tp*) and Ho Is the Initial mixture mass-specific enthalpy. 

Newton-Raphson corrector equations with log variable corrections (for 
self-scaling of the widely-varying mole numbers) are given by 


(A2) 


E 


3F (m) 3F (m) 

- LB Sr 4,09 3139 ’ - 1 - » 


^ 3,09 <C.*i 

N 


alogT- 


E 


0f( m ) a p( m ) 

— LULL 41o9 „<*> , T.n>1 , T (m> _ ( m ) 

k»n+l , v T (m) y n+1 T,n+1 


^ 3109 “M,l 


alog 


(A3) 


( A4) 


n+1 


3F 


The partial derivatives In equation (A3) are given by 

J 


(m) 

1 .n+1 


alog o£ m) , 
a k,n+l 


6 a (m) ^ 

= iy> + p Z-w (a U ” a lj )(R J a kj ' R -j a kj 


). i = l.N 


(A5) 


j=l 


and 
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3F 


(m) 


8i °9 ci " p 3=1 


J 


- R ,(B , + E ,/RT* m ] - nl ‘ + 1 


In the above equations, Is the Kronecker symbol 

= 0 ; 1 * k 


( A6) 


6 1k 


( A7) 


= 1 ; 1 = k 


and n'j and nj 1 are the molecularitles of the forward and reverse reactions 
j respectively, 


•i ■ I •« 


N 


(A8) 


"j' ' & “ii 


In equations (A3, A5 and A6), the partial derivatives with respect to the 
reciprocal of the mean molar mass, o m , are assumed to be negligible In com- 
parison with the other terms. 

The partial derivatives In equation (A4) can be derived from equation (A2) 
and are given by 


3F 


(m) 

T.n+1 


3 '°9 


rt (m) tr T (m) 
" °k,n+l n k n+1 


3F 


(m) 


N 


T.n+1 (m)^ (rr 

3log T<™> = '• 

n+ ' 1=1 


,(m) 7 T (m) 

n+1 pi n+1 


(A9) 


(A10) 


where c ,(T ,) Is the constant-pressure molal-speclf 1c heat of species 1 at 

pi n+1 

temperature . Substitution of equations (A9) and (A10) Into equation (A6) 

results In the following log-variable correction equation for the temperature 

N N 

(m) r # T ( m >iAinn J m > + T (m) V' rt (m) c fT (rn) Uloa T (m) - -F (m) 

Zj c k,n+l h k (T n+l )Alog °k,n+l + T n+1 2-w °1.n+l c p1 (T n+l )Alog T n+1 " F T,n+l 

k=l 1=1 

(All) 
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Equations (A3) and (All) are solved simultaneously by LU-decomposItlon and 
back-substitution, and the resulting log variable corrections are used to up- 
date the current Iterate values of n+1 ^n+1 by the approximate 

equations 


and 


(m+1 ) 
1 ,n+l 


(m) 

1 ,n+l 


1 + Alog a 


(m) 

1 ,n+l 


1 = 1 ,N 


T (m+1) (m) [ 

n+1 = n+1 [ 


* 41 »9 T n?l_ 


(A12a) 


(A12b) 


To start this Iteration process, the predicted values denoted by , (1 = l.N) 

i t n+l ' 

and are obtained quite simply by setting them equal to the values at the 

current time step 


a 


( 0 ) 

1 ,n+l 


l.n* 


1 = 1 ,N 


(A13a) 


T 


( 0 ) 

n+1 



(A13b) 


Jacobl-Newton Iteration 


The form of equations (Al), (A3), (A5) and (A6) was chosen to ensure di- 
agonal dominance of the Jacobian matrix. If It Is further assumed that the 
off-diagonal elements can be neglected with respect to the diagonal elements, 
equations (A3) and (All) can be rewritten as 


and 


aF (m) 

Lillil Aina a (m) p (ro) . < 1N 

dlog 409 W ’ 


T (m) y (m) ~ / T (m)w loa T (m) F (m) 

'n+1 fy °1,n+l c p1 (l n^l )A,0g T n+1 3 " F T,n+l 
In equation (A14), the partial derivatives are given by 


( Al 4) 


( Al 5) 


3F 


(m) 

1 ,n+l 


S1 °9 " -1 


which Is approximated by 


Jm) 

°1 .n+1 - 

U.h + * 


J=1 


= 1 ,N 


( Al 6) 
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3F 


(m) 


.(m) 


1 « n+1 -1 \ / i n , 1 1 p \ 

* ‘ = u,h 2-r u r j a u R -r 


ai °9 


( A1 7) 


M 


which when combined with equation (14) gives 


-p(m) (m) 

ah 1.n+l q 1 .n+1 n (m) 
al0Q Jm) = U h i »n+l 
8109 q 1 ,n+l 1 


(A18) 


where dI 11 ^ , , the destruction rate of species 1, can be evaluated along with 
1 , n+ 1 


f)’"' without calculating the entire Jacobian matrix. With this 


(m) 

1 ,n+l 

simplification, equation (A14) can be solved explicitly for the Iterative 
corrections 

/(^ * <*1-) =' ■ 1 •" <ai9 > 

The temperature corrector equation (A15) can be solved for the log- 
temperature correction 


Alnn ^ c(^) 

Alog T n+ -| = - F T ,n + V 


T (m) y (m) ~ / T (m)\ 

n+1 fa °1 ,n+l C p1 y n+1 J 


(A20) 


As with NR Iteration, the Iterate values of oi, n +l and T n+1 are 
updated using the approximate equations (A12) 


(m+1) (m) [", ^ (m) 

>1,n.l ‘ * 4,09 ,n+lj ’ 


1 = 1 ,N 


and 


(m+1) 
n+1 ' 1 


(m) f 

n+! [ 


1 * 41 og T<”> 


:] 


(A12a) 


(A12b) 


.( 0 ) 


To start the Iteration process, the predicted values are obtained 

from equation (9). (with 6^ = Z^) 

.( 0 ) 


0 ] ' , = a. + hf 
1 , n+1 1,n 


'exp(Z 1 #p h)-l 

'• n L _ 0 J : 


1 . 1,N (A21 ) 

The predicted temperature Is obtained by a single Newton-Raphson Iteration. 
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(A22) 


r ( 0 ) 

n+1 


T n + 


Ho-z «S“1 ,k,<v 


1=1 lf 


(0) c (T ) 
n+l c pv n' 


Convergence Test 

For both NR and JN Iteration, the test for convergence of the Iterates 
o^nlj Is based on the values alog and Is given by 



(A23) 


where, e Is a user-supplied local relative error tolerance. The above test 
Is used only for species whose mole numbers are greater than 10~20 ; ^ <e ., 
the summation does not Include species with mole numbers less than or equal to 
10"20. In addition, mole numbers less than 10 _ 20 are se t equal to 10 - 20. 

If convergence Is not obtained after ITMAX Iterations, where ITMAX Is the user- 
supplied maximum number of corrector Iterations to be attempted, the steplength 
Is reduced as discussed In the section Steplength Control and the step retried. 
If convergence Is achieved In M Iterations (M < ITMAX), the step Is accepted 
as successful and the solution Is updated. 


0 M*1 1 - 1 - N 


T - T<"> 
n+1 n+1 


(A24) 
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APPENDIX B 


AL60R1 - 

ALG0R2 - 
ALG0R3 - 
BSOLV - 

CONVG - 

CREKEQ - 
CREKO - 

CREK1D - 

DCOMP - 
DERIVS - 

EQUIL - 

ERATIO - 
GLOBL - 

HCPG - 

MECHN - 

NAMLSI - 

TEMPRK - 


DESCRIPTION OF SUBPROGRAMS 

Filters initial conditions. Computes Initial stepslze, Iterates 
L-formulated Implicit Euler approximations to convergence. 

Nonstiff regime solver. Employs JN Iteration to solve equations. 

Stiff regime solver. Employs NR Iteration. 

Standard routine for forward and back substitution of a previously 
LU-decomposed matrix. 

Manages control and monitoring of both NR and JN Iteration; Indi- 
cates when convergence criteria have been satisfied. 

Manages calls to subroutine EQUIL. 

Initializing routine for elemental and thermochemical data. Reads 
and catalogs data In NASA format from data file. 

Main routine. Sets Initial tlmestep, manages control of solution 
until end of prescribed tlmestep, and returns solution to calling 
program. 

Performs standard LU-decomposItlon of a square matrix. 

Evaluates kinetic expressions, and on demand, elements of Jacobian 
matrix. 

Calculates adiabatic flame temperature and equilibrium species dis- 
tribution for a mixture of gases at prescribed pressure and enthalpy. 

Calculates fuel-air equivalence ratio of a mixture of gases. 

Calculates kinetic rates and contributions to Jacobian for specially 
-prescribed global kinetic rate expressions. 

Evaluates enthalpy, constant-pressure specific heat capacity and 
one-atmosphere specific molar Gibbs function of a mixture of gases. 

Initializing routine for reading and cataloging kinetic rate data 
and establishing reaction stoichiometry vectors. 

Initializing routine for reading and cataloging problem control 
parameters, debugging options, etc. 

Performs a single Newton Iteration to determine the temperature of a 
given mixture of gases. 
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APPENDIX C 


SAMPLE INPUT AND OUTPUT 
Sample Input 


ELEMENTS 


c 

12.01115 

4.0 

H 

1.00797 

1 . 

0 

15.9994 

-2.0 

N 

14.0067 

0.0 


THERMO 

CO J 9/65C 1.0 1.00 0.00 O.G 300.000 5000.000 

0 . 2984068 9E 01 0 ..14891387E-Q2-0 . 57899678E-06 0 . 10364576E-09-0 . 69353499E-14 
-0.14245227E 05 0.63479147E 01 0.37100916E 0 1-0 . 16 1 ? *64E-02 0 . 36 923584E-05 
“0 . 20 31 967 3E-08 0 . 23953344E-12-0 . 1435630 9E 05 0.29555340E 01 
C02 J 9/65C 1.0 2.00 0.00 O.G 300.000 5000.000 

0.44608040E 01 0 . 30981717E-02-0 . 12392566E-05 0 . 2274 1 323 E-0 9-0 . 15525948E-13 
“0 . 48 96 1 4 38 E 05-0 . 98635978E 00 0.24007788E 01 0 . 87 35 0 9 0 5 E-02-0 . 66 07 086 1 E“05 
0 .20021860E-08 0. 63274039 E- 15-0. 48377520E 05 0.96951447E 01 
H J 9/65H 1.00 0.00 0.00 O.G 300.000 5000.000 

0.25000000E 01 0.0 0.0 0.0 0.0 

0.25471625E 05-0 . 46 0 1 1758E 00 0.25000000E 01 0.0 0.0 

0 . 0 0.0 0 . 2547 1 6 25 E 05-0. 46011 7 58 E 00 

H2 J 3/6 1H 2.0 0.0 0.0 O.G 300.000 5000.000 


0. 3100188 3E 01 
-0.87738013E 03- 

0 .55210343E-08- 
M20 

0.27167616E 01 
-0.29 9 05820 E 05 
-0 .296374Q4E-08 
N 

0.2450 26 78E 01 

0. 561160 35E 05 
-0 .56475602E-10 
NO 

0 . 3188 9992E 01 

0 . 98283242E 04 
-0.611 3 9254E-08 
N2 

0 . 28 96 3 1 94 E 01 
-0.90586182E 03 
-0 . 63217520 E- 0 9* 
0 

0 . 25420580E 01- 

0.29230801E 05 
-0.16028432E-08 
OH 

0.29106417E 01 

0 . 3935381 1 E 04 

0.18713971E-09- 

02 

0 . 3621 9521 E 01 
-0.12019822E 04 
-0 .67635071E-08 


0 .51119458E-03 
0 . 1 96294 12E 01 
0.18122726E-11- 
J 3/6 1 H 2.0 
0 .29451370E-02- 
0.66305666E 01 
0 .80702101E-12- 
J 3/6 IN 1.00 
0. 10661458E-03- 
0 . 44487 572E 01 
0 .20999038E-13 
J 6/63N 1.0 

0 . 13382279E-02- 
0.67458115E 01 
0.1591 9 072E-1 1 
J 9/65N 2.0 

0.15154863E-02- 
0.61615143E 01 
0 . 22577253 E- 12- 
J 6/620 1.00 

0.27550603E-04- 
0.49203072E 01 
0.38906 964E-12 
J 3/660 l.H 
0.95931627E-03- 
0.54423428E 01 
0.22571089E-12 
J 9/650 2.0 

0.736 18256 E-0 3- 
0.36150 942E 01 
0.21555977E-11- 


0 . 526442 0 4 E 
0. 30574446 E; 
0 . 988 9 0 4 3 0 E 
1.00 0.00 
0 . 80224 368 E 
0.40701275E 
0 .30279719E 
0.00 0.00 
0.74653315E 
0.25030699E 
0 . 56 0 988 98 E 
1.00 0.00 
0 . 528 9 931 6 E 
0.40459509E 
0 . 97 453867 E 
0.0 0.0 
0 . 57 23527 5E 
0 .36748257E 
0. 10611 587E 
0.00 0.00 
0.310 28 029E 
0. 2946423 3E 
0.29147641E 
1.00 0.00 
0.19441700E 
0 . 3837 5 9 3 1 E 
0 .36412820E 
0.0 0.0 
0 . 1 965221 9E 
0.36255 98 OE 
0 . 1 047 5225E 


-07-0 . 3490996 4 E- 10 
01 0 .26765198E-02- 
0 3- 0 , 22 9 97 0 46 E U1 
O.G 300.000 5000 

-06 0. 10226681E-09- 
01-0. 11084 4 99 E- 02 
05-0 . 3227 0 0 38E 00 
O.G 300.000 5000 

-07 0 . 18796520E-10- 
01-0 .21800181E-04 
05 0.41675749E 01 

O.G 300.000 5000 

-06 0 . 95919314E-10- 
01-0 .34181783E-02 
04 0 . 2 9 97 4976 E 01 
O.G 300.000 5000 

-06 0.99807385E-10- 
01-0.12081496E-02 

04 0 . 2358 0 4 18 E 01 

O.G 300.000 5000 

-08 0.45510670E-11- 

01-0. 16381664 E- 02 

05 0.29639931E 01 

O.G 300.000 5000 

-06 0.13756646E-10 

01-0.10778855E-02 
04 0.49370009E 00 
O.G 300.000 5000 

-06 0.3620 1556E-10- 
01-0. 18782183 E- 02 
04 0 . 43 0527 6 9E 01 


0.36945341E-14 

0.58099149E-05 

.000 

0 . 484721 04E-14 
0 .41521180E-05 

.000 

0 . 1 0259837 E~14 
0.54205284E-07 

.000 

0 . 64847 928E-14 
0 . 7 981 9174E-05 

.000 

0. 65223536 E- 14 
0 . 2324 0 1 0 0 E-05 

.000 

0 . 43680494 E-15 
0.24210303 E-05 

.000 

0.14224542E-15 

0.96830354E-06 

.000 

0.28945623E-14 
0 . 7 0554543 E-05 


MECHANISM 

CO OH C02 H 11.49 0.0 0.596 CGS 1 
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H 

02 


0 

OH 

H2 

0 


H 

OH 

H20 

0 


OH 

OH 

H 

H20 


H2 

OH 

N 

02 


NO 

0 

N2 

0 


N 

NO 

NO 

M 


N 

0 

H 

H 

M 

H2 


0 

0 

M 

02 


H 

OH 

M 

H20 


H2 

02 


OH 

OH 



14.34 

0.0 

16.492 

CGS 


13.48 

0.0 

9.339 

CG5 


13.92 

0.0 

18.121 

CGS 


14.0 

0.0 

19.870 

CGS 


9.81 

1.0 

6.250 

CGS 


. 13.85 

0.0 

75.506 

CGS 

M 

20.60 

-1.5 

149.025 

CGS 

M 

18.0 

-1.0 

0.0 

CGS 

M 

18.14 

-1 . 0 

0.340 

CGS 

M 

23.88 

-2.6 

0.0 

CGS 


13.0 

0 , 0 

43.0 

CGS 


& INPUT EPS=1 .QE-02,ITMAX=10,TKIN=1Q00.0,PATM=1Q.0,LDEBUG=.F. , 


NDEBUG=1,DELT=2. 0,SECS=1 . E-3, 5T0P=1 . E-0 3 SEND 
REACT ANTS 


C 1. 

r— t 

o 

CO 

1.0 

M 

G 

STOICH , 

H 2. 


H2 

2.0 

M 

G 

PYRLIZED 

N 2. 


N2 

7.52 

M 

G 

CH4-AIR 

0 2. 


02 

1.5 

M 

G 

MIXTURE. 
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H O sfl W vl CMJ1 ^ W IS) 



ELEMENTS 

C 12.011149 4.000000 
H 1.007970 1.000000 
0 15.999399 -2.000000 
N 14.006700 0.000000 


THERMO 

CO J 9/65 C 1.0 1.00 0.00 0. G 300.000 5000.000 0 

0 . 29840689E 01 0 . 148 91 387 E- 02-0 . 57899678E-06 0 . 1 0 36 4 57 6 E- 0 9- 0 . 6 93534 6 5E- 1 4 0 

-0.14245227E 05 0.63479147E 01 Q.371 00916E 0 1 - 0 . 1 6 1 9 0 96 4 E- 0 2 0 . 36 923575E- 0 5 0 

-0.20319673E-08 0 .23953344E-12-0 . 14356309E 05 0.29555340E 01 0 


C02 J 9/65 C 1.0 2.00 0.00 0. G 300.000 5000.000 0 

0.44608040E 01 0 . 30 981717 E-02-0 . 12392566 E-05 0 . 22741 323E-0 9-0 . 15525948E-13 0 

-0.48961438E 05-0 . 98635978E 00 0.24007788E 01 0 . 87350 905E-02-0 . 66 07 0852E-05 0 

0.20021860E-08 0 . 63274039E-15-0 . 48377520E 05 0.96951447E 01 0 


H J 9/65 H 1.00 0.00 0.00 0. G 

0.25000000E 01 0.00000000 0.00000000 0.00000000 

0.25471625E 05-0 .4601 1758E 00 0.25000000E 01 0.00000000 

0.00000000 0.00000000 0 . 25471625E 05-0 . 46011758E 00 


300.000 5000.000 

0.00000000 0 

0.00000000 0 

0 


0 


H2 J 3/61 H 2.0 0.0 0.0 0. G 300.000 5000.000 0 

0.31001883E 01 0 . 51 1 1 9458E-0 3 0 . 526 4420 4 E-0 7 - 0 . 349 0 9964E-1 0 0 . 36945341E-14 0 

-0.877 38 013E 0 3-0 . 1 96 294 12E 01 0.30574446E 01 0 . 26 7 6 5 1 98 E-0 2- 0 . 58 0 9 91 4 0 E-05 0 

0. 55210343 E- 08-0. 18122726E-11-0.98890430E 0 3- 0 . 22 9 97 046 E 01 0 


H20 

0.27167616E 01 
-0.29905816E 05 
-0. 29637404 E- 08 


J 3/61 H 2.0 1.00 0.00 0. G 300.000 5000.000 

0.29451370E-02-0.80224368E-06 0. 10226 68 1 E-0 9- 0. 4847 2070 E-14 0 

0 . 66305666E 01 0.40701275E 0 1 -0 . 1 1 084499E-02 0 . 4 152 1 180 E-05 0 

0 .80702101E-12-0 . 30279719E 0 5- 0 . 3227 0 0 32 E 00 0 


0 


0 . 24502678E 01 
0.56116031E 05 
-0. 56475602 E- 10 


J 3/61 N 1.00 0.00 0.00 0. G 300.000 5000.000 

0 . 10661458E-Q3-G . 746 533 15E-07 0 . 187 96 52 0 E-l 0-0 . 1 025 98 37 E-14 0 

0.44487572E 01 0.25030699E 0 1-0 . 2 1 8 0 0 18 1 E-0 4 0 . 5420 5284 E-0 7 0 

0.20999038E-13 0.56098895E 05 0.41675749E 01 0 


0 


NO 

0.31889992E 01 
0 . 982832 0 3 E 04 
-0 . 61139254 E-08 


J 6/63 N 1.0 1.00 0.00 0. G 300.000 5000.000 

0.1338227 9 E-02-0 .528 9 931 6 E-0 6 0 . 9591 9300E-1 0-0 . 64847928E-14 0 

0.67458115E 01 0.40459509E 0 1 -0 . 3418 1 783E-02 0 . 7981 9165E-05 0 

0.1591 9072E-1 1 0.97453828E 04 0.29974976E 01 0 


0 


N2 J 9/65 N 2.0 0.0 0.0 0. G 300.000 5000.000 0 

0 . 28 96 31 94 E 01 0 . 1515486 3E-02-0 . 57235275E-06 0 . 9980 7385E-1 0-0 . 65223536 E-14 0 

-Q.90586182E 03 0.61615143E 01 0.36748257E 0 1 - 0 . 12 08 1 4 96 E-0 2 0 . 2324 0 1 00 E-05 0 

-0. 63217498 E- 09-0. 22577253E-12-Q.10611587E 04 0.23580418E 01 0 


0 J 6/62 0 1.00 0.00 0.00 0. G 300.000 5000.000 0 

0 . 2542 0 580 E 0 1-0 . 2755 0 6 0 3E-04-0 . 31 028 027 E-08 0 . 4551 067 0 E-11-0 . 43680494E-15 0 

0.29230801E 05 0.49203072E 01 0.29464283E 01-0 . 16381664E-02 0 . 2421 0294E-05 0 

-0.160 28 432 E-08 0 . 38 9 0 6 96 4 E- 1 2 0.29147641E 05 0.29639931E 01 0 

OH J 3/66 0 l.H 1.00 0.00 0. G 300.000 5000.000 0 


Sample output 



0.29106417E 01 0.95931627E-03-0.19441700E-06 0 . 137 56 64 6 E-l 0 0 . 14224542E-15 0 

0.39353809E 04 0.54423428E 01 0.38375931E 0 1-0 . 10778853E-02 0 . 96830354E-06 0 

0 . 18713971E-09-0.22571089E-12 0.36412820E 04 Q.49370009E 00 0 

02 J 9/65 0 2.0 0.0 0.0 0. G 300.000 5000.000 0 

0 . 36219521E 01 0. 73618256 E- 03-0. 19652214 E~ 06 0 . 3620 1556 E-10-0 . 28945621E-14 0 

“0 . 12Q19822E 04 0.36150 942E 01 0.36255980E 01-0 . 18782183E-02 0 . 70554543E-05 0 

“0. 67635071 E-08 0 . 2155596 9E-1 1-0 . 1 0475225E 04 0.43052769E 01 0 


MECHANISM 


1 . CO 

OH 


C02 

H 




11.490 

0.000 

0.596 

CGS 1 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

CO EF 

= 

13.545 

0 .000 

22.716 

6 .513E-02 

9.996E-01 

2. H 

02 


0 

OH 




14.340 

0.000 

16.492 

CGS 2 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

13.235 

0.000 

0.476 

2.083E-02 

9.262E-01 

3. H2 

0 


H 

OH 




13.480 

0.000 

9.339 

CGS 3 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

13.123 

0.000 

7.308 

2.608E-03 

1.000E 00 

4. H20 0 


OH 

OH 




13.920 

0.000 

18.121 

CGS 4 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

12.923 

0.000 

1.025 

1 . 380E-Q2 

9.922E-01 

5. H 

H20 


H2 

OH 



14.000 

0.000 

19.870 

CGS 5 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

- 

13.359 

0.000 

4.803 

1 . 634E-02 

9.995E-01 

6 . N 

02 


NO 

0 




9.810 

1 .000 

6.250 

CGS 6 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

12.800 

0.000 

41.408 

5.705E-02 

9.999E-01 

7. N2 

0 


N 

NO 



13.850 

0.000 

75.506 

CGS 7 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

13.196 

0.000 

0.281 

4 .718E-03 

9 . 881E-0 1 

8 . NO 

M 


N 

0 


M 


20.600 

-1.500 

149.025 

CGS 8 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

- 

11.420 

0.000 

-6 .109 

5. 933E-02 

9.959E-01 

9. H 

H 

M 

H2 



M 


18.000 

-1 .000 

0. 000 

CGS 9 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

17 . 949 

0.000 

101 . 227 

5 . 105E-02 

1.000E 00 

10 . 0 

0 

M 

02 



M 


18.140 

-1 .000 

0 . 340 

CGS 10 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

18.837 

0 .000 

115.552 

7 . 434E-02 

1.000E 00 

11 . H 

OH 

M 

H20 



M 


23.880 

-2.600 

0 .000 

CGS 11 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

z 

18.611 

0 .000 

111.213 

1.489E-01 

9.999E-01 

12. H2 

02 


OH 

OH 



13.000 

0.000 

43.000 

CGS 12 


CALCULATED 

REVERSE RATE 

DATA, 

STD DEV 

AND 

CORR 

COEF 

= 

11.538 

0 .000 

24.955 

1.825E-02 

l.OOOE 00 

1. CO OH 


— 

= 

C02 

H 






2. H 

02 


— 

= 

0 


OH 






3. H2 0 


— 

= 

H 


OH 






4. H20 0 


— 

= 

OH 


OH 






5. H 

H20 


— 

= 

H2 


OH 






6. N 

02 


— 

: 

NO 


0 






7. N2 0 


— 

= 

N 


NO 








8. 

HO 


M 

= 

N 

0 

PI 

9. 

H 

H 

M 


H2 


PI 

10. 

0 

0 

m 


02 


PI 

11. 

H 

OH 

M 

r 

H20 


PI 

12. 

H2 

02 


= 

OH 

OH 







KINETIC RATE DATA 

IN SI UNITS 





J MODR 


ID 


BX 

TEN 

TACT 


BX2 

TEN2 

1. 

1 

1 10 

2 3 


8.490 

0.000 

299.949 


10.545 

0.000 

2. 

1 

3 11 

9 10 


11.340 

0.000 

8299.945 


10.235 

0.000 

3 . 

1 

4 9 

3 10 


10.480 

0.000 

4700 . 047 


10.123 

0 .000 

4 . 

1 

5 9 

10 10 


10.920 

0.000 

9119.770 


9.923 

0.000 

5 . 

1 

3 5 

4 10 


11 . 000 

0.000 

9999.996 


10.359 

0.000 

6 . 

1 

6 11 

7 9 


6.810 

1.000 

3145.446 


9.800 

0.000 

7 . 

1 

8 9 

6 7 


10.850 

0.000 

37999.977 


10.196 

0.000 

8 . 

2 

7 0 

6 9 


17.600 

-1.500 

74999.938 


5.420 

0.000 

9 . 

3 

3 3 

4 0 


12.000 

- 1.000 

0.000 


14.949 

0 . 000 

1 0 . 

3 

9 9 

11 0 


12.140 

-1.000 

171.112 


15.837 

0 . 000 

11 . 

3 

3 10 

5 0 


17.880 

-2.600 

0 .000 


15.611 

0 . 000 

1 2 . 

1 

4 11 

10 10 


10.000 

0 .000 

21640.668 


8.538 

0.000 

4INPUT 











DELT = 2.0 











EPS= 0 . 9999998E-02 










ITMAX= 10 











TKIN= 1000.0 











P ATM- 10.0 











TINY= 0.10E-19 










SECS = 0 . 9999999E-03 » 

29*0. 10 El 0 








STOP = 0.9999999E-03 










LDEBUG= F 











NDEBUG= 1 











t END 











REACTANTS 











C 1.00000 

0 1 . 

.00000 

0 .00000 


0.00000 CO 

1.00000 M G 

1 



H 2.00000 

0 . 

.00000 

0.00000 


0.00000 H2 

2.00000 M G 

1 



N 2.00000 

0 , 

.00000 

0 .0 0 0 0 0 


0.00000 N2 

7.52000 M G 

1 



0 2.00000 

0 , 

.00000 

0.00000 


0.00000 02 

1.50000 M G 

1 



*** REACTANT STREAM 

1 KKK 









I SPECIES 


MOLECULAR 

WEIGHT MOLE 

NUMBERS 

MASS FRACTIONS 






( KGMOLE I )/( KG 

I) (KGMOLE I)/(KG X) 

(KG I)/(KG X) 




1 . CO 



2.801E 

01 

3.440E-03 

9.636E-02 




2. C02 



4.401E 

01 


0.000 

0.000 




3. H 



1.008E 

00 


0.000 

0.000 




4. H2 



2.016E 

00 

6.880E-03 

1.3L7E-02 




5. H20 



1 .802E 

01 


0 .000 

0.000 




6. N 



1.401E 

01 


0.000 

0.000 





TACT2 

11432.273 

239.505 

3677.744 

515.668 

2417.458 

20839.703 

141.567 

-3074.641 

50944.621 

58153.891 

55970.461 

12558.973 



7. NO 

8 . N2 

9. 0 

10. OH 

11. 02 


3.001E 01 
2.801E 01 
1.600E 01 
1.701E 01 
3.200E 01 


0.000 
2.587E-02 
0 .000 
0 .000 
5.160E-03 


0 .000 
7.247E-01 
0.000 
0 .000 
1 .651E-01 


CO 

o 


TEMPERATURE 
ENTHALPY = 
PRESSURE = 
DENSITY = 
MEAN MOL UT 


1.000E 03 
5.090E 05 
1.013c 06 
2 . 947 E 00 
2.418E 01 


DEG K 
JOULES/KG 
N/M*#2 
KG/M# * 3 
KG/KGMOLE 


HSUBO, ER, PA , SM , RHOP, TK, TAU, TIME , N5TEP = 

5.090E 05 1.000E 00 1.013E 06 3.674E-02 1.266E 00 2.6192E 03 
SPECIES NAMES 

CO C02 H H2 H20 N 

02 

SPECIES MOLE NUMBERS, S2CI) 

6.6051E-04 2 . 77 94 E- 0 3 3.8895E-05 2.4178E-04 6.4945E-03 

2 . 8757 E-04 


0.000 9.990E 09 0 

NO N2 O OH 

8.1164E-09 1 . 96 47 E~ 0 4 2.5770E-02 


SPECIES MOLE FRACTIONS 


1.7976E-02 7 . 5645E-02 1.0586E-03 6.5802E-03 1.7675E-01 2.2090E-07 

7.8265E-03 

H5UB0, ER, PA, SM, RHOP, TK, TAU, TIME, NSTEP = 

5.090E 05 1.000E 00 1.013E 06 3.666E-02 1.269E 00 2.6186E 03 7.064E-05 1.00 
SPECIES NAMES 


5.3472E-03 


0E-03 99 


7 . 0136E-01 


CO C02 H H2 

02 

SPECIES MOLE NUMBERS, S2CI) 

H20 

N 

NO 

N2 0 

OH 

6.1412E-04 2 . 8292E-0 3 3.9434E-05 

3 . 9741 E-04 

2.1601E-04 

6.3588E-03 

1 . 1887E-08 

5.7265E-05 

2 . 5840E-02 

SPECIES MOLE FRACTIONS 






1.6752E-02 7.7173E-02 1.0756E-03 

5 . 8921E-03 

1 . 7345E-01 

3.2425E-07 

1 . 5620E-03 

7.0484E-01 


1 .0840E-02 


2.5897E-05 2.4847E-04 


7.0481E-04 6.7624E-03 


3.2643E-05 2.7593E-04 


S.9041E-04 7.5267E-03 
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Figure 3. - Variation of CPU time with local relative error tolerance for test problem 1. 
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Figure 4. - Variation of CPU time with local relative error tolerance for test problem 2. 
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